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ABSTRACT 

As the surface adopted for geodetic computation is a mathematical surface which is different from the physical surface, the geoid 
adopted as a reference for the vertical coordinate system, the ellipsoidal heights obtained from GPS obseryation are transformed 
to practical heights known as orthometric heights. The transformation of the ellipsoidal heights to orthometric heights requires 
the knowledge of the geoid-ellipsoid separation at the point of obseiyation. Since the geometric method requires the computation 
of geoid heights of points from GPS observation and geodetic leveling carried out over long distances which are labor intensive 
and prone to human errors, the accurate geoid heights of the points should be obtained from gravity measurement and a 
geometric geoid surface fitted to the gravimetric geoid heights. This paper presents detailed procedures for determining local 
gravimetric-geometric geoid model of an area or a region. The detailed procedures which consist of selection of suitable/evenly 
distributed points, DGPS and gravity observations of selected points, processing ofDGPS and gravity observations, computation 
of gravimetric geoid heights of the points, fitting of geometric geoid surface to the computed gravimetric geoid heights and 
computation of accuracy of the geoid model are presented in sequential order. 

Keywords: Gravimetric-geometric, Geometric geoid surface, Gravimetric geoid height, Stokes integral. 


1. INTRODUCTION 

The advent of Global Positioning System, GPS which determined accurately the three dimensional coordinates of points on the 
earth surface has led to the determination of local geoid models in various parts of the world. This is because the third dimension 
which is the ellipsoidal heights is determined on a mathematical surface as well as the ellipsoid. The ellipsoid is different from the 
surface adopted as a reference for vertical coordinate system. The surface adopted as a reference for vertical coordinate system is 
known as the geoid. The geoid is the surface which coincides with that surface to which the oceans would conform over the entire 
earth, if free to adjust to the combined effects of the earth's mass attraction (gravitation) and the centrifugal force of the Earth's 
rotation. Specifically, it is an equipotential surface, meaning that it is a surface on which the gravitational potential energy has the 
same value everywhere with respect to gravity (Olaleye et al, 2013). The transformation of the GPS ellipsoidal height to more 
meaningful or practical height known as the orthometric height requires the knowledge of the geoid-ellipsoid separation as well as 
the geoid height at the point of observation. The relationship among the ellipsoidal, orthometric and geoid heights according to 
Komarov (2007) is 

N = h — H (1) 

where, N is the geoid height, h is the ellipsoidal height and PI is the orthometric height. Figure 1 also shows the relationship 
among the ellipsoidal, orthometric and geoid heights with respect to the reference surfaces. 
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Figure 1: Relationship among Orthometric, Geoid and Ellipsoidal Heights. 

Source: Isioye and Youngu (2009) 

With the height difference determined using the orthometric heights, suitable gradient can be chosen and water can flow from one 
end to the other. 

The geoid can be determined using various methods such as the gravimetric, geometric, astro-geodetic and transformation 
methods amongst others. 

The gravimetric method can be carried out by the well-known Stokes-integral, equation (2) and the use of accurately determined 
absolute gravity data (Heiskanen and Moritz, 1967, Eteje, 2015 and Eteje et al, 2018). 

W = —ffAgS(^)r/<x (2) 

4 Try a 

Where, N is geoid undulation, A g is gravity anomaly, S{lf/ ) is stokes function, y is normal gravity on the reference ellipsoid 
and R is mean radius of the earth. The geometric method is to use the known “geoid heights” at some points, which are derived 
from co-located GPS-determined heights and levelled heights to interpolate the geoid heights at other points (Chen, and Luo, 
2004). The interpolation of the geoid heights at any other point involves the use of interpolation models such as bicubic model and 
other models like multiquadratic model, etc. In astro-geodetic method of geoid determination, the geoid heights of points are 
determined with reference to a geodetic (reference) station whose geoid height is known. The geoid heights differences between 
points are determined using the components of deflection of the gravity vector which can be obtained by carrying out 
astronomical and geodetic observations. The astronomical observation is carried out to determine the astronomical coordinates 
(astronomical latitude, <4* and astronomical longitude, A) by observing stars. The geodetic observation is used to determine the 
geodetic latitude, (p , geodetic longitude, A and ellipsoidal heights, h as well as the azimuths, OC and geodetic distances, t between 
network points. Using the astronomical and geodetic coordinates, the components of deflection of the gravity vector can be 
computed. The transformation method involves the use of the well-known Euclidean similarity transformation model which is 
used to convert Cartesian coordinates between two geodetic reference frames that generally differ in terms of three translation 
parameters (tx, ty, tz), three orientation parameters (ex, sy, ez) and a factor of uniform spatial scale change (5s). 

Using the geometric method, the geoid heights of points are obtained from GPS observations and geodetic levelling. The levelling 
of these points is carried out with respect to the Mean Sea Level, MSL as well as the geoid over long distances as these points are 
selected to cover the entire study area. Ideally, these points are selected along different roads so that their orthometric heights can 
be determined using closed loops differential levelling. Based on the fact that the points are chosen along various roads for 
accessibility and intervisibility reasons during levelling, these points are not evenly distributed over the entire study area. Also, the 
geodetic/differential levelling is labour intensive and prone to human errors over long distances. With the gravimetric method, the 
points can be chosen such that they are evenly distributed to cover the entire study area or region. This is because gravity 
observation does not require line of sight which in turn require intervisibility of points. Carrying out gravity observation for geoid 
heights computation is not labour intensive. Although, it is expensive but the utmost concern here is the accuracy with which the 
geoid heights are determined. Using the grid as well as the summation method of the modified Stokes integral in which the study 
area is divided into different compartments, the gravimetric geoid heights of each of the compartments are computed and 
averaged. The average geoid height is the geoid model of the study area. Taking the mean geoid height as the determined geoid 
model, it is certain that there will be residuals if the differences between the mean geoid and the geoid heights of the points are 
computed. This is because geoid heights vary with no constant value. 

Therefore, to determine accurate local geoid model of an area, the gravimetric and the geometric methods should be combined 
such that the geoid heights of the selected points are obtained from gravity observation and a geometric interpolation surface is 
fitted to the gravimetric geoid heights. Fitting an interpolation surface to the gravimetric geoid heights enables geoid heights of 
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new points to be obtained with the determined geoid model using variables such as geographic and rectangular coordinates. The 
fitting of the geometric geoid surface to the gravimetric geoid heights requires the model parameters to be computed. The 
computation is carried out with least squares adjustment method. 

This paper presents detailed procedures for determining local gravimetric-geometric geoid model of an area or a region. 

2. GRAVIMETRIC GEOID HEIGHTS COMPUTATION USING THE MODIFIED STOKES 
INTEGRAL 

Using the modified Stokes integral given in equation (1), the geoid height of points can be computed if their gravity anomalies and 
geographic coordinates are known. Featherstone and Olliver (1997) gave the integration of equation (1) as well as the Stokes 
formula as 
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(3) 


Where, N is the geoidal height of individual point, \j/ o is the surface spherical radius, y is the theoretical as well as normal 
gravity, A g is the gravity anomaly and r — R is the mean radius of the earth. 


2.1 Computation of Surface Spherical Radius, y/ 0 


The surface spherical radius, y/ n is computed as (Shrivastava et al, 2015) 

cos iff = sin $>sin (p x + cos (pc os (p' cosl/l 1 - A) (4) 

Where, 

(p = Mean latitude of the points 
(p l = Latitude of individual point 
A = Mean longitude of the points 
A 1 = Longitude of individual point 

2.2 Computation of Normal Gravity, y 


The normal as well as the latitude gravity is computed on a specified ellipsoid. According to Hinze et al (2005), the normal gravity 
is computed using: 


g e (\ + k sin 2 0) 
(l-e 2 sin 2 (j. i) 1/2 


(5) 


2 

Where, e~ is the first numerical squared eccentricity, g ( , is normal gravity at the equator and A: is a derived constant = 

0.001931851353. Based on equation (5), Aziz et al (2010) gave the formula for the computation of the normal or theoretical 
gravity as: 


S T 


9.7803267714 


' 1 + 0.00193185138639sin 2 (j) " 
V (1 - 0.00669437999013 sin 2 0) 1/2 J 


( 6 ) 


2.3 Computation of Gravity Anomaly, Ag 

The gravity anomaly is the difference between the points gravity reduced to the geoid and the latitude as well as the normal 
gravity computed on a specified ellipsoid and corrected for free air and effect of rock. The points gravity anomalies are grouped 
into two: the free air gravity anomalies and the bourgue gravity anomalies. The free air gravity anomaly, Ag j a is computed as 
(Aziz et al, 2010): 

^8 fa = g0bs~ S fa~ So ( 7 ) 
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Where, A g j a is free-air anomaly, g obs is observed gravity and g a is theoretical gravity and g j a = 0.3086[mGal/mJH. 

The bourgue gravity anomaly, A g B is given by Vermeer (2016) as 

Ag B =Ag fa -0.1119// (8) 

The orthometric heights of the points can be obtained from the Digital Terrain Model of the area. 

2.4 Computation of Mean Radius of the Earth, r = R 

The mean radius of the earth, R is computed using: 

R = - JMN (9) 

Where, M is the radius of curvature in meridian section and N is the radius of curvature in prime vertical.The formula for 
computation of the radius of curvature in prime vertical, /V is given as (Olaleye et al, 2013) 


N = ■ 


a 


(1 - e 1 sin 2 (p) ] 


GO) 


while that for computation of the radius of curvature in meridian section, M is given as (Kotsakis, 2008) 


M = ■ 


a( 1 -e~) 


(1- e 2 sin 2 (pf 


( 11 ) 


Where, 

a = semi-major axis 

(p = latitude of observation point 


e — 2/ — / 2 = eccentricity squared 
a-b 

J — -= flattening 

a 


b = semi-minor axis 


2.5 Computation of Combined Topographic Effect 


To obtain a precise geoid height of a point, the combined topographic effect is calculated and applied to the computed geoid 
height of the point. The formula for the computation of the combined topographic effect, SN^° b is given as (Sjoberg, 2000 and 
Kuczynska-Siehien et al, 2016): 


SN Topo - - 

UI v Comb 


2nGp 
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H 2 + — H 2 
3 R 


( 12 ) 


Were, G is the earth gravitational constant, p is density, R is the mean radius of the earth and H is the orthometric height of 
observation point which can be obtained from the DTM of the area. 


3. GEOMETRIC GEOID SURFACES 


Geometric geoid surfaces are mathematical interpolation surfaces fitted to geoid heights to enable geoid heights of new points to 
be determined using variable such as geographic or rectangular coordinates of the points. These surfaces include: plane surface, 
bi-linear surface, second degree surface, third degree polynomial and fifth degree polynomial. The surface to be adopted as well as 
the degree and order of the polynomial depends on the size of the study area and the variation of the geoid heights. For small area, 
the plane surface is used, for relatively large area, bi-linear surface is applied while for large area either the third or the fifth order 
polynomial surface is used. 

3.1 Plane Surface 

According to Alevizakou and Lambrou (2011), when the area of interest is small and the geoid has a normal variation, then it is 
possible to approach it by using a plane with a mean slope. The equation of such plan as given by Alevizakou and Lambrou 
(2011) is 

h, - //,■ = N fj (<p, ,A, ) = a 0 + a x .{(p i ~(p o ) + a 2 .(/l,. - ) (13) 

where, (p i , /l ( = geodetic coordinates, latitude and longitude, of point i. 

(p o ,X a - geodetic coordinates, latitude and longitude, of a central point in the region of interest. 
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N j = undulation of the geoid in the point i. 

h -, H t = geometric and orthometric heights in point i 

a o ,a 1 ,a 2 - unknown parameters of the plane. 

3.2 Bi-Linear and Second Degree Surfaces 

The bi-linear and the second degree surfaces are used where the study area show a more prominent surface. Alevizakou and 
Lambrou (2011) respectively gave the bi-linear and the second degree surfaces as 

^ - Hi = N fi (,Pi ,*,) = a 0 + a v {(Pi -p a ) + a 2 .(A, -A 0 ) + a 3 .(<Pi ~ PMA - A) 

^ - H t - N fi (Pi ,A) = a 0 + ip, ~p o ) + a 2 .(A, -A 0 ) + a 3 .(p, - p o f + 

+ a 4 .(A -Aj 2 +a 5 •(P , - p o )XA - A ) 


(14) 

(15) 


3.3 Polynomial Surface 

The polynomial surface used when determining geoid is given as (Kirici and Sisman, 2017) 

m k 

Ar u,r)=Z H a iJ xi y j (16 > 

i=0 j—k—i 
1=0 

Where, 

a j = polynomial coefficients 

m = degree of polynomial 
X, y = plan coordinates of point 

In applying the polynomial, the degree should be chosen and the polynomial should be formed for the chosen degree.Kirici and 
Sisman (2017) gave the 3rd degree polynomial as 

N = a ao + X + a 01 Y + $ 20 -^ + XY + cIq 2 Y ~ + ci 3 g X -\r Qt^X ~Y + r?j 2 XY~ + Qq 3 Y (17) 

Where, 

Y = ABS(y-y 0 ) 

X = ABS(x-x 0 ) 

y = Northing coordinate of observed station 
X = Easting coordinate of observed station 

y = Northing coordinate of the origin (average of the northing coordinates) 

X = Easting coordinate of the origin (average of the easting coordinates) 

The fifth degree polynomial which was applied in an area of 50 X 45 km~ in Turkey by Erol and (('elik (2004) is given as 

N(X,Y) = A () + \ n X + A W Y + \ 2 X 2 +A n XY + A 20 Y 2 +Aq 3 X 3 + A l2 X 2 Y + 

A 2l XY 2 + A 30 Y 3 + AuX 4 + A l3 X*Y + A 22 X 2 Y 2 + A 31 XY 3 + A 4Q Y 4 + (18) 

4, A 5 + A u X 4 Y + A 23 X 3 Y 2 + A 32 X 2 Y 3 + A 41 XY 4 + A 5Q Y 5 

Where, 

X = kx(p-p o ) 

Y = kx(y l- A) 

k = 100/p° ) = scaling and unit adjustment factor 
p — latitude of observation point 
A = longitude of observation point 
P a = latitude of the origin (average of the latitudes) 

A a = longitude of the origin (average of the longitudes) 
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3.4 Observation Equation Method of Least Squares Adjustment 


The fitting of geometric geoid surface to a set of geoid heights required the model parameters to be computed. The computation of 
these parameters is done using observation equation method of least squares adjustment. The functional relationship between 
adjusted observations and the adjusted parameters as given by Ono et al (2014) is: 

L a = F(X a ) (19) 


Where, L a = adjusted observations and X a = adjusted parameters. Equation (19) is linear function and the general observation 

equation model was obtained. The system of observation equations is presented by matrix notation as (Mishima and Endo, 2002 
and Ono et al, 2018): 

V = AX-L (20) 

But the residual matrix is not necessary when applying least squares adjustment technique for determination of local geometric 
geoid model parameters. Thus the general matrix notation becomes 

L = AX ( 21) 

where, A = Design Matrix, X = Vector of Unknowns, L = Observation Matrix. 

X={A t A)-'A t L (22) 

The step by step procedures for determining geometric geoid model parameters are detailed in Eteje and Oduyebo (2018). 


3.5 Evaluation of the surface fitting 


In choosing a geometric geoid surface to be fitted to a particular region, the fitting of the surface has to be evaluated to determine 
if the chosen surface is actually fitting to the region. According to Alevizakou and Lambrou (2011), the unknown parameters must 
be tested in order to assess whether they are statistically important for a certain confidence level (for example 95%). Each of the 
unknown parameters a t must conform to the following formula (Alevizakou and Lambrou, 2011): 

"^95% — a i ^3) 

Where, <7 a = standard deviation of each unknown parameter Cl i 

Z 95% = coefficient of the normal distribution for one dimension arrays for confidence level 95%. 

The parameters a t will be considered statistically important only if equation (23) is valid. The standard deviation of the unknown 
parameters, <7 is obtained from the variance-covariance matrix generated during least squares computation of the model 
parameters. 


3.6 Accuracy/Reliability of Gravimetric-Geometric Geoid Model 


The accuracy of the determined local gravimetric-geometric geoid model is obtained using the Root Mean Square Error, RMSE 
index. To evaluate the determined local gravimetric-geometric geoid model accuracy, the local geoid model is used to determine 
the geoidal heights of points whose gravimetric geoid heights are known. The gravimetric-geometric geoid model geoidal 
undulations are compared with the known gravimetric geoidal undulations of the points to obtain the residuals. The Root Mean 
Square Error, RMSE index for the computation of gravimetric-geometric geoid model accuracy as given by Kao et al (2017) and 
Eteje and Oduyebo (2018) is 

\v T v 

RMSE = ±J - (24) 

V n 


Where, 


v = N ggh - N g _ gmgh (Residual) 


N GGH = Gravimetric Geoid Height of Point 
Xg-gmgh ~ Gravimetric-Geometric Model Geoid Height of Point 
n = Number of Points 
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4. PROCEDURES FOR LOCAL GRAVIMETRIC-GEOMETRIC GEOID MODEL 
DETERMINATION 

The local gravimetric-geometric geoid model determination involves the following procedures: selection of suitable/evenly 
distributed points, DGPS and gravity observations of selected points, processing of DGPS and gravity observations, computation 
of gravimetric geoid heights of the points, fitting of geometric geoid surface to the computed gravimetric geoid heights and 
computation of accuracy of the geoid model. Figure 1 shows the local gravimetric-geometric geoid model determination 
procedures flow chart. 



Figure 2: Local Gravimetric-Geometric Geoid Model Determination Procedures Flow Chart 

4.1 Selection of Suitable/Evenly Distributed Points 

Select suitable points within the study area such that the points are evenly distributed to cover the entire area. This will enable the 
gravimetric geoid heights of the study area to be evenly determined to cover the whole area. 

4.2 DGPS and Gravity Observations of Selected Points 

Carry out DGPS and gravity observations to respectively obtain the geographic coordinates and the gravity of the selected points. 
The geographic coordinates of the points are used during the processing of the gravity observations. The gravity observation of the 
points should be carried out with a gravimeter. 

4.3 Processing of DGPS and Gravity Observations 

The acquired DGPS and gravity observations should be processed on the local ellipsoid adopted for geodetic computation in the 
study area. The normal gravity, y should be computed on the local ellipsoid using equation (5). The free air and bourgue gravity 

anomalies, A g of the points should be respectively computed using equations (7) and (8). 
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4.4 Computation of Gravimetric Geoid Heights of the Points 

Having computed the gravity anomalies of the points, the next step is to compute the gravimetric geoid heights of the points using 
equation (3). To apply equation (3) for computation of the geoid heights of the points, the mean radius of the earth, r and the 
spherical radius, \j/ 0 have to be respectively computed using equations (9) and (4). The computation of the mean radius of the 

earth requires the computation of the radius of curvature in prime vertical, N and radius of curvature in meridian section, M using 
equations (10) and (11) respectively. To obtain precise geoid heights of the points, the combined topographic effect should be 
calculated with equation (12) and applied to the computed geoid heights of the points. 

4.5 Fitting of Geometric Geoid Surface to the Computed Gravimetric Geoid Heights 

Having computed the precise gravimetric geoid heights of the points, a geometric geoid surface is fitted to the geoid heights to 
enable geoid heights of new points to be determined within the study area using their rectangular or geographic coordinates. Any 
of the geometric geoid surfaces given in equations (13) to (18) can be used depends on the size of the study area and the variation 
of the gravimetric geoid heights of the points. It is advisable to fit two or more geometric geoid surfaces to the geoid heights so as 
to determine the surface that is most suitable for application in the study area. To fit a geometric geoid surface to the gravimetric 
geoid heights, the model parameters, CL have to be computed using least squares technique as well as equation (22). To assess 
whether each of the computed model parameters is statistically important for a certain confidence level, equation (23) is applied. 

4.6 Computation of Accuracy of the Geoid Model 

The determined gravimetric-geometric geoid model has to be validated to determine its reliability. The determination of the 
reliability of the geoid model requires the computation of the accuracy of the model using the Root Mean Square Error, RMSE 
index as well as equation (24). To do this, the differences between the gravimetric geoid heights of randomly selected points for 
validation and their corresponding model geoid heights are computed to obtain the residual matrix, V in equation (24). 

5 CONCLUSION 

The accurate transformation of geometric heights obtained from GPS observation to practical as well as orthometric heights has 
resulted in the determination of local geoid models in different parts of the world. Local geoid models are determined using 
various methods such as gravimetric, geometric, astro-geodetic, astro-gravimetric, transformation and hybrid methods. The 
geometric method requires the computation of geoid heights of selectedpoints from GPS observation and geodetic as well as 
differential levelling carried out over long distances which is labour intensive and prone to human errors.Therefore, the accurate 
geoid heights of thepointsshould be obtained from gravity measurement and a geometric geoid surface fitted to the gravimetric 
geoid heights. Consequently, the this paper has in simple terms presented detailed procedures of fitting a geometric geoid surface 
to gravimetric geoid heights. Thus, local gravimetric-geometric geoid model determination. 
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